/*******************************************************************************

                                APODUTIL

  Title:        Utilities for running apodfc and apodfo. 

  Author:       Hanna Szoke

  Date:         10/5/00: full comparison of 2 wil files for scaling 
		fobs to fcalc; see scale_fo().

*******************************************************************************/

#include "util.h"
#include "cellparams.h"
#include "dims.h"	
#include "apodize.h"

	/* variables declared in apodize.h	*/

int	Nbin ;		/* # bins used in program */
real	delta ;		/* smearing factor */
real	binwidth ;	/* width of a bin */
real	*w_lnIsq ;	/* wilson y coords */
real	*apodwt ;	/* weights associated with Iave */

	/* variables local to apodutil.c */

static	real	*Ibin ;		/* Binned, not averaged, I's */
static	real	*sigsqbin ;	/* Binned, not averaged, sigma's */
static	real	*Iave ;		/* Binned & averaged I's */
static	real	*sigsqave ;	/* Binned & averaged sigma's */
static	int	*cases ;	/* # of entries in a bin */

static	int	lin_orig ;	/* index of start of linearized signal */
static	int	m0 ;		/* index of first non-empty bin */
static	real	min_x ;		/* min x corresponding to max data res */
static	real	max_x ;		/* max x corresponding to min data res */
static	real	y_intercept ;
static	real	slope ;

static	real	min_res = 0.0 ;	/* linearization range minimum */
static	real	max_res = 0.0 ;	/*       "        "    maximum */
static	real	apod_res = 0 ;

static	void	send_it_out_asc(char *, real *, int, float, float) ;
static	void	generateBins(int, real *, real *, real *) ;
static	char	checkBins() ;
static	void	averageBins() ;
static	void	newvar() ;
static	int	read_it_in_asc(char *, real **, real **) ; 
static	int	genlinBins(real *) ;		

#define	BINWID 0.002 	/* default width of a bin */

void a_readInput()    /* Read input, identifying data in terms of keywords. 
                         Do requisite reciprocal lattice calculations.  */
{
     int	k ;

     readBasicInput() ;
     binwidth = BINWID ;

     for (k = 0; k < Nlines; k++) {
          strcpy(nextline, allinp + k*MAXSTRING) ;
          sscanf(nextline, "%s", id) ;

          if (strcmp(id, "APOD_RES") == 0) {
               sscanf(nextline, "%*s %g", &t1) ;
               apod_res = t1 ;
          }
          else if (strcmp(id, "BINWIDTH") == 0) {
               sscanf(nextline, "%*s %g", &t1) ;
               binwidth = t1 ;
          }
          else if (strcmp(id, "MIN_RES") == 0) {
               sscanf(nextline, "%*s %g", &t1) ;
               if (t1 == 0)
                    Fatal("Illegal min_res") ;
               min_res = t1 ;
               min_x = 1.0 / (min_res*min_res) ;
          }
          else if (strcmp(id, "MAX_RES") == 0) {
               sscanf(nextline, "%*s %g", &t1) ;
               if (t1 == 0)
                    Fatal("Illegal max_res") ;
               max_res = t1 ;
               max_x = 1.0 / (max_res*max_res) ;
          }
          else if (strcmp(id, "NUSQ_RANGE") == 0) {
               sscanf(nextline, "%*s %g %g", &t1, &t2) ;
               if ((t1 <= 0) || (t2 <= 0))
                   Fatal("Illegal Nusq range") ;
               min_x = t1 ;
               max_x = t2 ;
               min_res = 1./ sqrt(t1) ;
               max_res = 1./ sqrt(t2) ;
          }
     }

     if (apod_res == 0)
	  apod_res = input_res;

     if (min_res < max_res) {
         Warning( "Your resolution limits seem to be inverted!") ;
	 Fatal(" MAX_RES should be a smaller number than MIN_RES") ;
     }

     if ((min_res == 0) ||(max_res == 0)) {
        printTwice("Input values for Nu squared range are missing.") ;
        printTwice("Code will use data limits for linearization.") ;
     }
     else {
        sprintf(message, "Linearization limits are (%g, %g) %s", 
	     min_res, max_res, units) ;
        printTwice(message) ;
     } 

     sprintf(message, "Apodization resolution = %g %s", apod_res, units) ;
     printTwice(message) ;

     return ;
}

void	init_apod_arrays(int Ntot, real *Int, real *ksq, real *sigsqI)
{
     int	n ;

     for (n = 0; n < Ntot; n++) {
          *(Int+n) = 0 ;
          *(ksq+n) = 0 ;
          *(sigsqI+n) = 0 ;
     }
     return ;
}
void	allocate_Bin_arrays(int N) 
{
     cases   = (int *) e_malloc(N*sizeof(int), caller) ;
     Ibin = (real *) e_malloc(N*sizeof(real), caller) ;
     sigsqbin = (real *) e_malloc(N*sizeof(real), caller) ;
     Iave = (real *) e_malloc(N*sizeof(real), caller) ;
     sigsqave = (real *) e_malloc(N*sizeof(real), caller) ;
     w_lnIsq = (real *) e_malloc(N*sizeof(real), caller) ;
     apodwt = (real *) e_malloc(N*sizeof(real), caller) ;
}
void	do_binning(int Ntot, real *Int, real *ksq, real *sigsqI)
{
     char okay ;

     do {
          generateBins(Ntot, Int, ksq, sigsqI) ;
          okay = checkBins() ;
     } while (!okay) ;

     averageBins() ;
     newvar() ;

     return ;

}
			/*********************************************
              		Move Int, ksq, sigsqI from lists to bins.
			*********************************************/

void	generateBins(int Ntot, real *Int, real *ksq, real *sigsqI)
{
     int	m, n ;
     int	overflow = 0 ;
     real	Iover = 0, sover = 0 ;

     /*  Prepare bins */

     for (m = 0; m < Nbin; m++) {
	  *(Ibin+m) = 0 ;
	  *(sigsqbin+m) = 0 ;
	  *(cases+m) = 0 ;
     }

     for (n = 0; n < Ntot; n++) {

	  m = (*(ksq+n) + 0.5*binwidth) / binwidth ;

          if (*(Int+n) > 0) {

               if (m >= Nbin) {
	            overflow++ ;
	            m = Nbin-1 ;
                    Iover += *(Int+n) ;
	            sover += *(sigsqI+n) ; 
	       }

	       *(Ibin+m) += *(Int+n) ;
	       *(sigsqbin+m) += *(sigsqI+n) ;
	       *(cases+m) += 1 ;
	  }
     }

    if (overflow > 0) {
	  sprintf(message, 
	  "%d entries overflowed the highest resolution bin.", overflow) ;
	  printTwice(message) ;
	  sprintf(message, 
  "Their contribution to that bin was I = %g (sig = %g)", Iover, sover) ;
	  printTwice(message) ;
	  sprintf(message, 
  "                            out of I = %g (sig = %g)", *(Ibin+Nbin-1), 
							*(sigsqbin+Nbin-1)) ;
	  printTwice(message) ;
     }

     return ;
}

			/*********************************************
              		Check the binning procedure 
			*********************************************/

char	checkBins()
{
     int	m1, mlast=0 ;
     char	okay = TRUE ;

     /****************************************************
     Check bins: empty ones at the beginning are okay;
     empty ones at the end are also okay, but they cause 
     Nbin to be suitably redefined.   
     NOTE: an empty bin can occur either because no (hkl)
     entries go into it or because the total intensity,
     Iave (not yet averaged) is 0.
     *****************************************************/

     if (verbose)  {
	   fprintf(fp_log, "Initial bin information:\n\n") ;

           for (m0 = 0; m0 < Nbin; m0++) 
               fprintf(fp_log," m = %d, Ibin = %g, cases = %d\n",
	       m0, *(Ibin+m0), *(cases+m0)) ;

               fprintf(fp_log, "\n") ;
          }

     for (m0 = 0; m0 < Nbin; m0++) 
          if ((*(cases+m0) > 0) && (*(Ibin+m0) > 0))
              break ;

     for (mlast = Nbin-1; mlast >= m0; mlast--) 
          if ((*(cases+mlast) > 0) && (*(Ibin+mlast) > 0))
              break ;
		
     for (m1 = m0; m1 < mlast; m1++) 
          if ((*(cases+m1) == 0) || (*(Ibin+m1) == 0)){
               sprintf(message,
	       "There is at least one empty bin (# %d) within range (%d, %d)!",
	       m1, m0, mlast) ;
	       printTwice(message) ;
               okay = FALSE ;
               break;
          }	       

     if (!okay) {
          binwidth *= 2;
          Nbin /= 2.;
          sprintf(message, "Halving the number of bins to %d.", Nbin) ;
          printTwice(message) ;
     }
     else {
          fprintf(fp_log, "First and last occupied bins are %d and %d\n", 
	  m0, mlast) ;
          Nbin = mlast + 1 ;
     }
     return(okay) ;
}
			/*********************************************
              		Average bin contents based on # of cases.
			Note sigsqave is averaged with 1/sq(# of cases)
			*********************************************/
void	averageBins()
{
     real	d ;
     int	m ;
     char	bad = FALSE ;

     for (m = m0; m < Nbin; m++) {
	  *(Iave+m) = 0. ;
	  *(sigsqave+m) = 0. ;
     }
     
     for (m = m0; m < Nbin; m++) {
	  d = (real) *(cases+m) ;
	  *(Iave+m) = *(Ibin+m) / d ;
	  *(sigsqave+m) = *(sigsqbin+m) / (d*d) ;
     }
     
     if (verbose) {
	  fprintf(fp_log, "\n") ;
	  sprintf(message,
	     " bin      d(%2s)      1/d*d    count    ave(Isq)\n", units) ;
	  fprintf(fp_log, "%s", message) ;
	  fprintf(fp_log, " ---     -----     -------    -----    --------\n") ;

          for (m = m0; m < Nbin; m++)  {
	       d = 1. / sqrt((m+0.5)*binwidth) ;
	       fprintf(fp_log, "%4d  %8.2f  %10.5f  %5d  %10.2f\n", 
			    m,  d, (m+0.5)*binwidth, *(cases+m), *(Iave+m)) ;
               }
     }

     for (m = m0; m < Nbin; m++) {
	  if (*(sigsqave+m)  <= 0) {
	       sprintf(message, "sigsq array is not > 0 at m = %d", m) ;
	       Warning(message) ;
	       bad = TRUE ;
          }
     }
     if (bad)
	  Fatal("Quitting ...") ;
     
     return ;
}
void	newvar() 
{
     int	m ;

     for (m = m0; m < Nbin; m++) {
          *(w_lnIsq+m) = log(*(Iave+m)) ;
	  *(apodwt+m) = *(Iave+m) * *(Iave+m) / *(sigsqave+m) ;
     if (verbose)
	  fprintf(fp_log, 
	  "m = %d, ln Isq = %g, weight = %g\n", m, *(w_lnIsq+m), *(apodwt+m)) ;
     }
     return ;
}
int	genlinBins(real *lin_lnIsq)
{
     real	x ;
     int	m, m1 ;

     for (m = 0; m < Nbin; m++) {
	  x = (m+0.5)*binwidth ;
	  *(lin_lnIsq+m) = y_intercept + x*slope ;
     } 

     for (m1 = 0; m1 < Nbin; m1++) {
	  x = m1*binwidth ;
	  if (x > max_x) 
	       break ;
     } 

     /* set lin_orig */

     for (lin_orig = 0; lin_orig < Nbin; lin_orig++) 
	  if ((x = lin_orig*binwidth) > min_x)
	       break ;
     
     return(m1-lin_orig) ;
}

     /************************************************************************
     Derive best fit of Isq to y = y0 + slope*x, in range (min_x, max_x).
     See Numerical Recipes, 1986, p.505 (14.2.6).
     *************************************************************************/

real     linfit(int N, real *lnIsq, real *awt) 
{
     real	Sx = 0, Sy = 0, Sw = 0 ;
     real	Sxx = 0, Sxy = 0 ;
     real	Den = 0 ;
     real	x, y, w ;
     real	ksqmin = 100., ksqmax = 0 ;
     real	Bcryst ;
     int	m ;

     /*  Ensure non-zero value of max_x (if not set by user */

     if (max_x == 0)
          max_x = 1./dsq_limit ;

     for (m = 0; m < N; m++) {

	  x = (m+0.5)*binwidth ;
	  y = *(lnIsq+m) ;
	  w = *(awt+m) ;

	  if ((x > min_x) && (x < max_x)) {
	       Sw += w ;
	       Sx += x*w ;
	       Sy += y*w ;
	       Sxx += x*x*w ;
	       Sxy += x*y*w ;
          }
	  if (x < ksqmin)
	       ksqmin = x ;
	  if (x > ksqmax)
	       ksqmax = x ;
     }
     Den = Sw*Sxx - Sx*Sx ;

     if (Den == 0) {
          sprintf(message, "ksq range is (%g, %g), ", ksqmin, ksqmax) ;
          Warning(message) ;
          sprintf(message, 
            "min_x = %g (min_res = %g), max_x = %g (max_res = %g).",
            min_x, min_res, max_x, max_res) ;
          Warning(message) ;   
          Warning(
 "You should examine your values for min_res, max_res and/or binwidth\n") ;
          Fatal("Quitting...") ;
     }

     y_intercept = (Sxx * Sy - Sx * Sxy) / Den ;
     slope       = (Sw * Sxy - Sx * Sy) / Den ;

     fprintf(fp_log, 
       "\nCurve fitted over (%g, %g) to (y = y0 + slope*x)\n",  min_x, max_x) ;
     fprintf(fp_log, "with y0 = %g, (variance: %g)\n", y_intercept, Sxx / Den) ;

     fprintf(fp_log, "and slope = %g (variance: %g).  Cov(y0, b) = %g.\n", 
		slope, Sw / Den, -Sx / Den) ;

     /************************************************************************
     We use the slope of ln(Isq) vs. ksq to calculate Bcryst; the definition
     of Bcryst is:
	Bcryst = [2*pi*delfac*input_res]^2 * eta,
     where delfac and eta are defined in basics.c and are
	0.7 and 0.6 for BCC, and
	0.6 and 0.8 for SC.
     *************************************************************************/

     Bcryst = -2. * slope ;
     sprintf(message,
	"\nThe average crystallographic B factor is  %g %s sq,", Bcryst, units) ;
     printTwice(message) ;

     sprintf(message,
	"corresponding to a resolution of %g %s.\n", 
	     1./(TWOPI*grid_spacing/input_res) * sqrt(Bcryst/eta), units) ;
     printTwice(message) ;

     return (Bcryst) ;
}

real	get_delta(real Bcryst)
{
     real	ratsq ;

     /************************************************************************
     We use Bcryst to calculate delta, (where ksq stands for h*h/a*a + ...):
     It is found by solving:

   exp[-Bcryst*ksq/2]*exp[-2PISQ*delta*eta*drsq*ksq]=
                      exp[-2PISQ*eta*drsq*(apod_res/input_res)^2*ksq]  
     
     -Bcryst/2 - 2PISQ*delta*eta*drsq = -2PISQ*eta*drsq*(apod_res/input_res)^2,
     or
     delta = (apod_res/input_res)^2 - Bcryst / (4*PISQ*eta*drsq).

     In these formulas, drsq = dr*dr, dr stands for the resolution of the 
     Gaussians in our representation (grid_spacing).

     *************************************************************************/

     ratsq = (apod_res / input_res) * (apod_res / input_res) ;

     delta = ratsq - Bcryst / (4 * PISQ * eta * grid_spacing * grid_spacing) ;

     if (delta > 0) 
          sprintf(message, 
       "The data should be smeared using an apodizing factor of %g.\n", delta) ;
     else { 
          sprintf(message,  
       "The data should be unsmeared (apodizing factor is %g).\n", delta) ;
	  delta = 0 ;
     }
     printTwice(message) ;

     return (delta) ;
}

void make_wilson(real *dstdev)  
{
     char	out_filename[MAXSTRING] ; /* name of output file for x/y data */
     real	*l_lnIsq ;
     real	*d_lnIsq ;
     real	x, dave = 0 ;
     float	origin, interval ;	/* signal origin & interval */
     int	Nlin ;			/* length of lin_lnIsq */
     int	count = 0 ;
     int	m ;

  /*******************************************************
  Now prepare for display the original binned log(I^2),
  the linearly fitted version and the difference file.
  *******************************************************/
     
     l_lnIsq = (real *) e_malloc(Nbin*sizeof(real), "apodutil") ;
     d_lnIsq = (real *) e_malloc(Nbin*sizeof(real), "apodutil") ;

     Nlin = genlinBins(l_lnIsq) ;

     interval = binwidth ;
     origin = m0*interval ;

     sprintf(out_filename, "wil") ;
     send_it_out_asc(out_filename, w_lnIsq+m0, Nbin-m0, origin, interval) ;
     origin = lin_orig*interval ;

     sprintf(out_filename, "lin_wil") ;

     send_it_out_asc(out_filename, l_lnIsq+lin_orig, Nlin, origin, interval) ;

     for (m = 0; m < Nbin; m++)
          *(d_lnIsq+m) = *(l_lnIsq+m) - *(w_lnIsq+m) ;

     for (m = 0; m < Nbin; m++) {

	  x = (m+0.5)*binwidth ;

	  if ((x > min_x) && (x < max_x)) {
               dave += *(d_lnIsq+m) ;
	       count++ ;
          }
     }
     dave = dave / count ;
     *dstdev = 0 ;

     for (m = 0; m < Nbin; m++) {

	  x = (m+0.5)*binwidth ;

	  if ((x > min_x) && (x < max_x)) 
               *dstdev += (*(d_lnIsq+m)  - dave) * (*(d_lnIsq+m)  - dave) ;

     }
     *dstdev = sqrt(*dstdev / count) ;

     sprintf(message, 
     "Standard deviation of wilson curve with respect") ;
     printTwice(message) ;
     sprintf(message, "   to its linearized version is %g.", *dstdev) ;
     printTwice(message) ;

     e_free(Nbin*sizeof(real), l_lnIsq) ;
     e_free(Nbin*sizeof(real), d_lnIsq) ;

     return ;
}

void make_selected_wilson(char *name)  
{
     float	origin, interval ;	/* signal origin & interval */

     interval = binwidth ;
     origin = m0*interval ;

     send_it_out_asc(name, w_lnIsq+m0, Nbin-m0, origin, interval) ;

     return ;
}
void show_wilson_plots() 
{
     if (graphics) {
          printTwice("\n") ;

	  strcpy(message, 
	  "xmgrace wil lin_wil -legend load -timestamp") ;
          system(message) ;
     }
     else {
	  printTwice(
"You should use a graphics program such as xmgrace to plot wil and lin_wil") ;
     }
}

		/*******************************************************
 		Write an ASCII file of x/y info.
		*******************************************************/

void send_it_out_asc(char *out_name, real *array, int N, float x0, float delx)
{
     FILE 	*fp ;
     int	i ;

     if ((fp = fopen(out_name,"w")) == NULL) 
	  Fatal("Couldn't open ascii file for writing.") ;

     fprintf(fp, "@ title \"Wilson Plot\"\n") ;
     fprintf(fp, "@ xaxis label \"Nu_sq (1/res_sq) (%s^-2)\"\n", units) ;
     fprintf(fp, "@ yaxis label \"ln(F_sq)\"\n") ;

     for (i = 0; i < N; i++)		/* write pairs of x and y */
          fprintf(fp,"%g\t%g\n", x0+i*delx, *(array+i)) ;
	  
     fclose(fp) ;
}
		/*******************************************************
 		Read an ASCII file of x/y info. Return actual length, 
		Ntotal, which may be <= N, but not > N!
		*******************************************************/

int	read_it_in_asc(char *in_name, real **xarray, real **yarray) 
{
     FILE 	*fp ;
     int	i ;
     int	Ntotal = 0 ;
     real	*xap, *yap;

     check_exist(in_name) ;
     fp = fopen(in_name,"r") ;

     while (fgets(nextline, MAXSTRING, fp) != NULL)  
	  Ntotal++ ;
	  
     fclose(fp) ;

     *xarray = (real *) e_malloc(Ntotal*sizeof(real), caller) ;
     *yarray = (real *) e_malloc(Ntotal*sizeof(real), caller) ;
     xap = *xarray ;
     yap = *yarray ;

     for (i = 0; i < Ntotal; i++, xap++, yap++) {
	  *xap = 0 ;
	  *yap = 0 ;
     }
     fp = fopen(in_name,"r") ;
     i = 0 ;
     xap = *xarray ;
     yap = *yarray ;

     while (fgets(nextline, MAXSTRING, fp) != NULL)  {
          sscanf(nextline, "%s", id) ;
          if (strcmp(id, "@") != 0) {
             sscanf(nextline, "%g %g", &t1, &t2) ;
	     *(xap+i) = t1 ; 
	     *(yap+i) = t2 ;
	     i++ ;
          }
     }
	  
     fclose(fp) ;
     return(i) ;
}
real	scale_fo(char *name_fo, char *name_fc) 
{
     real	*xarr_fc, *yarr_fc ;
     real	*xarr_fo, *yarr_fo ;
     real	fcwidth, fowidth ;
     real	rat_oc, epsi ;
     real	diff = 0 ;
     real	correction ;
     float	origin, interval ;
     char	in_filename[MAXSTRING] ;
     int	Nfc, Nfo ;
     int	nc_skip = 1, no_skip = 1 ;
     int	nc_start = 0, no_start = 0;
     int	n, nc,no;
     int	count = 0 ;
     char	good ;

     sprintf(in_filename, "%s", name_fc) ;
     Nfc = read_it_in_asc(in_filename, &xarr_fc, &yarr_fc) ;
     sprintf(in_filename, "%s", name_fo) ;
     Nfo = read_it_in_asc(in_filename, &xarr_fo, &yarr_fo) ;

     /*******************************************************
     Scaling fits two _wil plots which may have different 
     origins, lengths, binwidths ... these must be checked.
     First: do their ranges overlap?
     Then: are their binwidths compatible?
     *******************************************************/

     if ((*xarr_fo >= *(xarr_fc+Nfc-1)) ||
         (*xarr_fc >= *(xarr_fo+Nfo-1))) {
          sprintf(message, "file1 limits: %g %g, file2 limits: %g,%g",
	          *xarr_fo, *(xarr_fo+Nfo-1), *xarr_fc, *(xarr_fc+Nfc-1)) ;
          printTwice(message) ;
          Fatal("Wilson files do not have overlapping ranges!") ;
     }

     fcwidth = (*(xarr_fc+1) - *xarr_fc) ;
     fowidth = (*(xarr_fo+1) - *xarr_fo) ;
     rat_oc = fowidth /fcwidth ; 
     epsi = 0.01 * fcwidth ;

     if (rat_oc >= 1.) {
         for (n = 1, good = FALSE; n <= 16; n++)
             if (fabs(fowidth - n*fcwidth) < epsi) { /*to avoid roundoff with floats*/
                 good = TRUE ;
                 nc_skip = n ;
                 break ;			/* widths are compatible */
             }
     }
     else {
         for (n = 1, good = FALSE; n <= 16; n++)
             if (fabs(fcwidth - n*fowidth) < epsi) { /*to avoid roundoff*/
                 good = TRUE ;
                 no_skip = n ;
                 break ;			/* widths are compatible */
             }
     }

     if (!good) {
	  sprintf(message, "fo binwidth is %g, fc binwidth is %g", 
		  fowidth, fcwidth) ;
	  printTwice(message) ;
	  Fatal("Fo and fc Wilson arrays are incompatible") ;
     }

     /* find first common point in xarr_fc and xarr_fo */

     for (nc_start = 0; nc_start < Nfc; nc_start++) {
	 for (no_start = 0; no_start < Nfo; no_start++) {
	     if ((*(xarr_fo + no_start) >= *(xarr_fc + nc_start))) 
	          break ; 
          }

	  if (fabs((*(xarr_fo + no_start) - *(xarr_fc + nc_start))) < epsi)
	       break ;
     }
     if ((no_start > Nfo/2) || (nc_start > Nfc/2))
	   Fatal(
	   "Couldn't find suitable common starting index for scaling.") ;

     for (nc = nc_start, no = no_start; (nc < Nfc && no < Nfo); 
          nc+= nc_skip,  no+= no_skip) {

          diff += *(yarr_fc + nc) - *(yarr_fo + no) ;
	  count++ ;
     }

     correction = diff / (float) count ;
     origin = *xarr_fo ;
     interval = fowidth ;

     sprintf(message, "\tProposed value for fscale is %g", 
             exp(correction / 2.)) ;
     printTwice(message) ;

     for (n = 0; n < Nfo; n++) 
	  *(yarr_fo + n) += correction ;

     send_it_out_asc("scaled_wil", yarr_fo, Nfo, origin, interval) ;

     if (graphics) {
        sprintf(message, "xmgrace %s %s -legend load -timestamp", 
                name_fo, name_fc) ;
        system(message) ;

        sprintf(message, "xmgrace %s %s -legend load -timestamp", 
                "scaled_wil", name_fc) ;
        system(message) ;
     }
     else {
	sprintf(message, "You can use a graphics program to plot %s and %s",
		"scaled_wil", name_fc) ;
        printTwice(message) ;
     }
     return(correction/2.) ;
}
